rm(list = ls())
require(mfx)
require(zTree)
require(stargazer)
require(stats4)
require(lmtest)
require(bbmle)
require("multiwayvcov")
require(msm)
require(nlWaldTest)
require(sandwich)
require(miceadds)


zTT1= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October20ChangingPriors.xls")
#zTT1= zTreeTables("C:/Users/nneligh/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October20ChangingPriors.xls")
subjects1=zTT1$subjects

zTT2= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October23ChangingPriors.xls")
#zTT2= zTreeTables("C:/Users/nneligh/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October23ChangingPriors.xls")
subjects2=zTT2$subjects
subjects2$Subject=subjects2$Subject+100

zTT3= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October24ChangingPriors.xls")
#zTT3= zTreeTables("C:/Users/nneligh/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October23ChangingPriors.xls")
subjects3=zTT3$subjects
subjects3$Subject=subjects3$Subject+200

subjects=rbind(subjects1,subjects2,subjects3)

priors=c(0.5,0.85)

#Clean data
subjects$responder=(subjects$"response1[50]">0)
subjects=subset(subjects,responder>0)

uidlist=unique(subjects$Subject)

#find N
#**
length(uidlist)
#**


#create a frame with observations
obsframe=data.frame(uid=vector('numeric'),block=vector('numeric'), priorA=vector('numeric'), number=vector('numeric'), state=vector('numeric'), response=vector('numeric'))
for(i in 1:length(uidlist)){
playframe=subjects[i,]
uidplace=uidlist[i]

#block 1
for (j in 1:50){
blockplace=1
priorplace=priors[playframe[1,682]]
numberplace=j
stateplace=playframe[1,(735+j)]
responseplace=playframe[1,(71+j)]
addframe=data.frame(uid=uidplace,block=blockplace, priorA=priorplace, number=numberplace, state=stateplace, response=responseplace)
obsframe=rbind(obsframe,addframe)
}

#block 2
for (j in 1:50){i
blockplace=2
priorplace=priors[playframe[1,683]]
numberplace=j
stateplace=playframe[1,(785+j)]
responseplace=playframe[1,(221+j)]
addframe=data.frame(uid=uidplace,block=blockplace, priorA=priorplace, number=numberplace, state=stateplace, response=responseplace)
obsframe=rbind(obsframe,addframe)

}

}

expframe=obsframe

#definitions
expframe$priorB=1-expframe$priorA
expframe$Aresponse=(expframe$response==1)
expframe$Bresponse=(expframe$response==2)
expframe$stateA=(expframe$state==1)
expframe$stateB=(expframe$state==2)
expframe$correct=(expframe$response==expframe$state)

write.table(expframe, "C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Revision/Reports/AmbhuelChangingPriors.csv", sep = ",", col.names = T, append = F)

#Begin constructing first table
pAgiven2=vector('numeric')
restrictionpAgiven1=vector('numeric')
pAgiven1=vector('numeric')
prob=vector('numeric')

#Base Regression
modelCP=lm(Aresponse~factor(stateA+priorA)+0,data=expframe)


#Adjust COV
vcovmodCP=vcovHC(modelCP, type="HC3")
vcovmodCPcluster=cluster.vcov(modelCP, expframe$uid)

modelcpsterr=vector("numeric")
#Model CP sterrors
for(i in 1:4){
modelcpsterr[i]=sqrt(vcovmodCPcluster[i,i])
}

logitCP=glm.cluster(data=expframe, Aresponse~factor(stateA+priorA)+0, cluster="uid", weights=NULL, subset=NULL,  family="binomial" )

prob=vector('numeric')
sterr=vector('numeric')
tstat=vector('numeric')

pAgiven2[1]=modelCP$coefficients[1]
pAgiven2[2]=modelCP$coefficients[2]

pAgiven1[1]=modelCP$coefficients[3]
pAgiven1[2]=modelCP$coefficients[4]

restrictionpAgiven1[1]=((2*priors[1]-1)/priors[1])+pAgiven2[1]*((1-priors[1])/priors[1])
restrictionpAgiven1[2]=((2*priors[2]-1)/priors[2])+pAgiven2[2]*((1-priors[2])/priors[2])

sterr[1]=as.numeric(deltamethod(~x3-((2*0.5-1)/0.5)+x1*((1-0.5)/0.5),modelCP$coef,vcovmodCPcluster))
tstat[1]=(pAgiven1[1]-restrictionpAgiven1[1])/sterr[1]
prob[1]=1-pnorm(tstat[1])

sterr[2]=as.numeric(deltamethod(~x4-((2*0.6-1)/0.6)+x2*((1-0.6)/0.6),modelCP$coef,vcovmodCPcluster))
tstat[2]=(pAgiven1[2]-restrictionpAgiven1[2])/sterr[2]
prob[2]=1-pnorm(tstat[2])

NIASframe=data.frame(priorA=priors,pagiven2=pAgiven2,restriction=restrictionpAgiven1,pagiven1=pAgiven1,prob=prob)

#*****************************
stargazer(NIASframe,rownames=FALSE,summary=FALSE)
#*****************************
